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SUMMARY 


A fourth-order box method is presented for calculating numerical solutions 
to parabolic, partial-differential equations in two variables or ordinary dif- 
ferential equations. The method, which is the natural extension of the second- 
order box scheme to fourth order, is demonstrated with application to the incom- 
pressible, laminar and turbulent, boundary-layer equations. The efficiency of 
the present method is compared with other two-point and three-point higher order 
methods, namely, the Keller box scheme with Richardson extrapolation, the method 
of deferred corrections, a three-point spline method, and a modified finite- 
element method. For equivalent accuracy, numerical results show the present 
method to be more efficient than the other higher order methods for both laminar 
and turbulent flows. 


INTRODUCTION 

Some of the recent effort in the application of numerical methods to the 
boundary-layer equations has been directed toward finite-difference methods 
which have truncation errors that are of higher order than the second-order 
methods in use today. The term higher order refers to the truncation error in 
the coordinate normal to the body surface. The truncation error in the tangen- 
tial coordinate is second order in all the methods discussed in this report. 

The advantages of higher order methods are twofold. First, they can be used to 
obtain a numerical solution as accurate as a second-order method with less com- 
puter time and storage; alternatively, they can be used to produce a signifi- 
cantly more accurate solution with the same amount of computer time and storage 
as a second-order method. 

The higher order methods can be either composite or intrinsic. Composite 
methods rely on one or more second-order calculations to achieve their higher 
order accuracy. Richardson extrapolation and the method of deferred correc- 
tions, both discussed by Keller (ref. 1), when applied to second-order methods 
are examples of composite methods. Intrinsic methods solve a single set of 
governing difference equations which have a truncation error smaller than 
second order and whose solution does not depend on a previous second-order 
calculation. 

The higher order intrinsic methods proposed for the boundary-layer equa- 
tions consist of finite-difference schemes which involve two or three grid 
points. The three-point schemes fall into two classes. The first class con- 
sists of procedures which are fourth order for uniform grids. These procedures 
treat the functional value and certain derivatives as unknowns at three colloca- 
tion points and can be derived via Taylor series or polynomial interpolation. 

In this category are the Fade approximation of Kreiss (ref. 2), or the so- 
called compact scheme (ref. 3), the Mehrstellen formulation (ref. 4), and the 
formulation of Peters (ref. 5). The second class of three-point schemes con- 
sists of methods for variable grids such as the spline collocation methods of 



Rubin and Graves (refs, 6 and 7) and Rubin and Khosla (refs. 8 and 9) • Adam 
(refs. 10 and 11) and Kordulla (private communication) have extended the com- 
pact and Mehrstellen schemes, respectively, to variable grids. The two classes 
of three-point schemes are similar in that the resulting finite-difference equa- 
tions involve three nodal points but are different in that the first class is 
restricted to uniform grids whereas the second class is applicable to variable 
grids. The emphasis in this study is on methods applicable to variable grids. 

One disadvantage of higher order intrinsic methods involving three nodal 
points is that the usual boundary conditions for incompressible flow (i.e., 
no slip and no injection at the surface and the tangential velocity component 
approaches the edge value as the boundary layer merges with the mainstream) are 
no longer sufficient, since the resulting general system of finite-difference 
equations contains more unknowns than equations. This difficulty is usually 
finessed by inventing one or more additional boundary conditions at the outer 
edge of the layer and by applying an additional equation at the surface bound- 
ary. The choice of which additional boundary conditions to use is not clear. 

Another disadvantage of some higher order intrinsic methods is the com- 
plexity of the resulting nonlinear, finite-difference equations and the associ- 
ated difficulty in solving them efficiently at each column. For example, the 
spline 4 method of Rubin and Khosla (refs. 8 and 9) would seem to require the 
solution of a 5 X 5 block- tridiagonal matrix to solve the fully coupled, incom- 
pressible, boundary-layer equations. A simpler solution scheme, which lags the 
continuity equation one iteration behind the momentum equation, is reported by 
Rubin and Khosla (ref. 8 ). Since the equations in this scheme are solved uncou- 
pled, the errors diminish in a linear manner with each iteration. In contrast, 
the better second-order methods, such as the Davis coupled scheme (ref. 12) and 
the Keller box scheme (ref. 13), solve the equations fully coupled with Newton 
linearization, and thus, for laminar flows, quadratic convergence is achieved. 
Hence, the advantages of some higher order methods, relative to second-order 
methods, may be reduced or lost entirely in practical engineering calculations 
because of slow convergence of the nonlinear system. 

The main purpose of this report is to present a fourth-order intrinsic box 
scheme for obtaining numerical solutions to parabolic, partial-differential 
equations or ordinary differential equations, with application here to the 
steady, two-dimensional, incompressible, laminar or turbulent, boundary-layer 
equations. The method has the following features: (1) It results in finite- 

difference equations that involve only two nodal points and therefore is for- 
mally fourth-order accurate on all grids, ( 2 ) it results in a 3 x 3 matrix of 
unknowns at each nodal point when the equations are solved in a coupled manner, 
(3) it utilizes Newton linearization and demonstrates quadratic convergence for 
laminar flows, and (4) it requires only the standard incompressible boundary 
conditions given previously. In short, the method is the natural fourth-order 
extension of the second-order box scheme. It is an example of a general class 
of high- accuracy, two-point methods examined by White (ref. 14) and discussed 
briefly by Keller (ref. 15). Keller offers an operation-count analysis that 
suggests that such methods may be superior to the well-known Keller box scheme 
(with Richardson extrapolation) to achieve acpuracy less than or equal to six 
order. 
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The second purpose of this report is to compare the present fourth-order 
box scheme with other proposed or existing higher order methods, including a 
second-order method with Richardson extrapolation, the method of deferred cor- 
rections, a three-point spline scheme, and a modified finite-element scheme. 


AnfBn 

^ij 

C 

C 


SYMBOLS 

3x3 matrices defined by equation (27) 
i,jth element of Aj^ 
i,jth element of 

constant in equation (37) ; also 3x3 matrix given by ecjuation (B2) 
constant in equation (38) 


Cf surface skin-friction coefficient 

^OO^X) j 

oq,c\ constants in grid-stretch function (see eq. (35)) 

Sl/C 2 fC 3 constants in finite-difference expression for 9( )/9^ terms 
(eq. (AD) 


^n 


3x3 matrix given by equation (B8) 


,(D .(2) .(3) 
■n r^n r ^n 

(D ^(2) ^(3) 


coefficients in finite-difference solution technique 
(see eqs. (B3) to (B5) ) 


absolute numerical error in the wall shear, percent 

Ej* absolute numerical error in boundary-layer displacement thickness, 

percent 

F damping factor (see eq. (10c)) 

f normalized boundary-layer velocity component in x-direction, u/U^ 

g represents any sufficiently differentiable quantity 

h step size in n “coordinate, equal to Arin-l 

I number of intervals across boundary layer 

k Von Karman constant, 0.41 

L* reference length 
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Z scale length measured in Y-coordinate 

I scale length measured in n -ccxjrdinate 

a = 1 + e 

N number of grid points across boundary layer 

components of vTOtor (see eq. (26)) 

Rp correction term (see eq. (33b) ) 

Rj, Reynolds number based on L*, P^U^L*A'<5 

Rx Reynolds number based on x*, P<jU*x*/Vi* 

s = 3f/3n 

T central-processor-unit time 

t defined by equation (10b) 

U* inviscid flow velocity in x-direction 

U = 

' 00 

u* time-averaged viscous flow component in x-direction 

u = 


friction velocity, 

V* time-averaged viscous flow component in Y-direction 

V = \IRlV*/j* 

V transformed viscous flow component in y-direction (see eq. (7)) 
vector defined by equation (26) 

X* position coordinate measured along body surface from leading edge 

or stagnation point 

X = x*/L* 

y stretched normal coordinate, 

y* physical distance normal to body surface, y = 0 at surface 

y = y*/L* 
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law-of-the-wall coordinate, u<Jy*/v* 

Zj^ vector defined by equation (25) 

bt exponent in grid-stretch function (see eq. (35c)) 

coefficients used in appendix A in finite-difference equations 
6 pressure-gradient parameter (see eq. (8) ) 

Ari step size in n -coordinate, = Tin “ ^n-1 

step size in ^-coordinate, A5j„_*l = “ ^m-1 

boundary-layer displacement thickness measured in n-coordinate 
(see eq. (34)) 

e nondimens ional eddy viscosity (see eq. (10a)) 

^n parameter in grid-stretch function defined by equation (35b) 

ri transformed normal coordinate (see eq. (5)) 

0 parameter determining where differential equation is evaluated in 

^-coordinate, 6=1/2 for Crank-Nicolson scheme and 0=1 for 
three-point backward difference scheme 

X mixing length (see eq. (lOd)) 

y* molecular viscosity 

V* kinematic viscosity, y*/P* 

5 transformed surface coordinate (see eq. (5)) 

p* density 

( 0u*\ 

T* wall shear, ly* j 

\ 3y*/y=0 

0 ) ccjefficient denoting truncation error (see eq. (37)) 

Subscripts: 

e inviscid flew conditions at y = 0 

exact value obtained with 640 intervals across boundary layer 
m grid index in ^-coordinate for finite-difference formulation 
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n 


grid index in n-coordinate for finite-difference formulation 


w viscous flow conditions at y = 0 

00 free-stream conditions 
Superscripts: 

1 present iteration 

’ differentiation with respect to ri-coordinate 

quantity evaluated outside main iteration loop 
* dimensional quantity 

-1 inverse of matrix 

vector 

T transpose of a matrix (or vector) 

Acronyms : 

B4S fourth-order box scheme 

CKBS conservative Keller box scheme 

CKBSRE conservative Keller box scheme with Richardson extrapolation 

DCS Davis coupled scheme 

KBS Keller box scheme 

MDC method of deferred corrections 

MFE modified finite-element method 

RKSl Rubin-Khosla S^(4,0) spline method 

FOURTH-ORDER BOX SCHEME (B4S) 

Governing Equations 

The present method is demonstrated with application to the steady, two- 
dimensional, incompressible, laminar or time-averaged turbulent, boundary-layer 
equations (ref* 16, p. 545) which are expressed here in Gortler variables 
(ref. 17): 
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Mcaientum 


LLH 

9n\ 3ti 


- V 


8f 

in 


- 3(f2 - 1) - 2Kf 


9f 

_ _ 0 


Continuity 


8v 9f 

— + 25 — + f 

3n 35 


0 


Boundary conditions are 

f=v=0 at n=0 (no slip, no injection) 
f -*■ 1 as n ~ 


Other quantities are given by 


5 = r Ue((|)) d(l) n 
^0 


Ug (x) 



u 

f = — 
Ue 


^ _ 2Cf 3n/3x 

V = V + 


U. 


U, 


e - — — 

Ue d5 

«, = !+£ 

The eddy viscosity e is given by reference ISs 
|3f| 


e = t 


3n 


where 


t = 


F2\2j|R^u; 

fc 


( 1 ) 


( 2 ) 


( 3 ) 

( 4 ) 


( 5 ) 


( 6 ) 


( 7 ) 


( 8 ) 

( 9 ) 


(10a) 


(10b) 
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p = 1 



(10c) 


X / k y\ 

- = 0.085 tanh (k = 0.41) (lOd) 

I \^0.085 ij 


Note that the variables in equations (1) to (10) are all nondimensional . 


Method of Solution 


The finite-difference expression around which the present method is formu- 
lated is given by 

h , , h^ „ „ 

9n - 9n-l “ -<9n + 9^-1) + ~ (9n ~ 9n-l) + 0(h5) =0 (11) 

where g represents any sufficiently differentiable quantity# h is the vari- 
able step size in the normal coordinate h (i.e.# b = - t)jj_i), and the 

primes denote differentiation with respect to ri. As mentioned previously, 
equation (11) is not new. Liniger and Willoughby (ref. 19) analyzed the sta- 
bility of equation (11) as an implicit method for solving initial-value prob- 
lems for stiff systems of ordinary differential equations. Equation (11) is the 
fourth-order member of a general class of high-accuracy, two-point equations 
discussed by White (ref. 14) and Keller (ref. 15). Hirsh (ref. 3) used the 
expression to formulate boundary conditions for the three-point compact scheme 
applied to the incompressible driven cavity problem. 


A key step in the present procedure, as in the second-order Keller box 
scheme, is to reformulate the problem in terms of a first-order set of partial- 
differential equations. This is done by defining s = 3f/9ri and rewriting 
equations (1) and (2): 


3 3f2 , 

— (Jls - vf) = 2K + (1 + 3)f2 - 3 

3n 3C 


( 12 ) 


3f 

= s (13) 

3n 


3v 

3h 



(14) 
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If a vector 


-► 

g 


is now defined as 


g 



(15) 

(16) 
(17) 


then g' is the right-hand sides of equations (12) to (14) and g" , their 
derivatives with respect to the n-coordinate; 


25 


9f2 

w 


+ (1 


+ 3)f2 - B 


g' 


- 


S 


> 


( 18 ) 

(19) 





( 20 ) 




r 


3 (fs) 

25 — - + (1 + B)fs 

35 


(21) 


J" = <(2Jl - D-l 


vs + B (f 2 - 1) +5 t'|s|s 


( 22 ) 


ds 

-25 s 

35 


(23) 






where equation (1) has been used to express s' in terms of the dependent vari- 
ables f, s, and V (eq. (22)). Substitution of equations (15) to (23) into 
equation (11) leads to three, nonlinear, finite-difference equations. The 
3( )/35 terms in these equations are written so that they may be approximated 
either by a two-point central-difference quotient (the Crank-Nicolson scheme) 
or by a three-point backward difference method; in either case, the equations 
are second-order accurate in the 5~coordinate and fourth-order accurate in the 
T) -coordinate. These nonlinear equations are then linearized by Newton's method. 
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with the exception of the term t'|s(s in equation (22) and the quantities I 
and Si which are used to compute X and y"*". Application of Newton's method 
to linearize the eddy viscosity terms was found to particularly accelerate con- 
vergence for the turbulent similarity cases. Newton linearization yields three, 
linear, finite-difference equations which can be expressed as 


^n^n-1 ®n^n ~ ^^n 

where 


^n ~ \®m,n»^m,n»%,n 
'^n “ 



an 

312 

313 ' 


’t>ll 

bi2 

bi3 

An - 

321 

322 

323 

Bn = 

t>21 

t>22 

t>23 


J31 

332 

333 


_t>31 

t>32 

t>33_ 


(24) 


(25) 

(26) 


(27) 


The superscript i in equation (25) denotes the present iteration value and 
the subscripts m and n are grid indices in the and r) -coordinate, respec- 
tively. The matrices and and the vector are functions of the 

dependent variables evaluated at the i-1 iteration and/or at the previous 
^-stations. Equations (24) must be solved repeatedly until an acceptable level 
of convergence is obtained. 

Although the double subscripted matrix elements ajj and b^j change 
values with the index n, it is suppressed for simplicity. Superscript i 
and subscript m will also be suppressed for the same reason, ^pendix A 
shows hew Af,, B^j, and are obtained. 

The boundary conditions are 

fl = V] = 0 (no slip, no injection) (28a) 

fN = T (28b) 

It is noteworthy that, since equations (24) can be applied at the outer 
boundary, the total number of finite-difference equations exactly balances the 
total number of unknowns. In contrast, a three-point method formulated in 
terms of f, s, and v could not be applied at the outer boundary and would 
result in more unknowns than equations; hence the method would require addi- 
tional boundary conditions and equations applied at the boundaries. 
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The solution of the resulting linear system of finite-difference equa- 
tions (24) has been given before; see Cebeci and Smith (ref. 20). The method 
used here is given in appendix B. 


HIGHER ORDER METHODS WITH WHICH B4S IS COMPARED 

A classification of the higher order methods examined in this study is 
given in table I. The methods in the table are described in the following 
sections. 


TABLE I.- CLASSIFICATION OF HIGHER ORDER METHODS 


Method 

Type 

Number of grid 
points in 
finite-difference 
equations 

Formal accuracy with - 

Uniform 

grids 

Variable 

grids 

B4S 

Intrinsic 

2 

Fourth order 

Fourth order 

RKS1 

Intrinsic 

3 

Fourth order 

Third order 

MFE 

Intrinsic 

3 

Fourth order 

Third order 

CKBSRE 

Composite 

2 

Fourth order 

Fourth order 

MDC 

Composite 

2 

Fourth order 

Fourth order 


Intrinsic Methods 

Three-point spline methods .- Two three-point spline methods were examined 
in this study. Both methods are formally fourth-order accurate for uniform 
grid spacing and can also be shown to be fourth order for analytical variable 
grid stretchings. They are the spline 4 method of Rubin and Khosla (ref. 8)^ 
which solves the momentum and continuity equations in an uncoupled manner , and 
the spline S^ (4^0) method, which solves the equations coupled. The equations 
for the S^ (4,0) method were given by Rubin and Khosla (ref. 9) without applica- 
tion. However, since this method is easily applied to solve the coupled momen- 
tum and continuity equations with dependent variables f , s = f • , and v and 
thus achieves quadratic convergence for laminar flows, the three-point spline 
method used for comparison here is the S^(4,0) method. Numerical results 
obtained with the spline 4 the (4,0) methods were virtually identical. 

Rubin-Khosla spline S^(4^0) (RKS1) .- The RKSl method is shown here applied 

to the momentum equation in similarity form which is obtained by deleting the 
3( )/9C term in equation (1). To apply the spline method to turbulent flow, 
the momentum equation is written in the following nonconservative form: 


fi,f" + (2.' - v)f ' - 3 (f2 - 1) = 0 


(29) 


n 



In the RKSl method, the f" derivative in the momentum equation is 
replaced by a spline relation involving the dependent variables s and f 
at three node points. The resulting finite-difference momentum equation is 
then coupled with finite-difference expressions for the continuity equation 
and an auxiliary equation which enforces spline continuity; both are obtained 
via spline relations. 

After linearizing the finite-difference momentum equation by Newton's 
method with the exceptions for eddy viscosity noted in the B4S, these 3x3 
block-tridiagonal equations are solved simultaneously for s^, fn» and v^ 
with boundary conditions, 

fl = VI = 0 (no slip, no injection) (30a) 

fN = 1 (30b) 

Since three-point finite-difference equations can be applied only at interior 
node points, there are three more unknowns than equations. However, vjj does 
not appear in the finite-difference equations written in nonconservative form, 
so that only two more unknowns than equations remain. The two additional equa- 
tions needed are obtained by setting sjj = 0 and evaluating s-j by applying 
equation (11 ) at n = 2 to the momentum equation (12). 

Modified finite-element method (MFE) .- The modified finite-element method 
of Rubin and Khosla (ref. 21) differs from the three-point spline methods in 
that the equations are solved in conservative form. The momentum finite- 
difference equation is given by 


9n+l “ 9n-l 


h 

2 


<^gA+l + (1 + Cf)gn + 9n-l 


* - (o2 - ,)g; - giUl] ■ 0 (31) 

where o = Arin/ATijj_i and g„, g^/ and gj^ are given by equations (15), 

(18), and (21). The continuity equation is the same as the B4S finite- 
difference expression (eq. (11)). These equations are coupled with an auxil- 
iary equation which enforces spline continuity. 


After linearizing the momentum equation with the exceptions for eddy vis- 
cosity noted in the B4S, these 3x3 block- tridiagonal equations are rewritten 
with dependent variables f , v, and M = f" subject to the boundary condi- 
tions given in equations (30) . 


The coupled system of equations for fjj» Vn» and Mjj written in this 
form is unstable when numerical solutions are attempted. However, a stable 
system can be obtained by eliminating and from the momentum dif- 

ference equation via the continuity equation. A computer code listing which 
supplied details of this method was provided courtesy of Rubin and Khosla. 
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since Vjj no longer appears in the difference equations, only two addi- 
tional equations are required. The momentum equation is used to show that 
M*| = - 3 - A spline relation relating Sj^ to % and fjj was also used 
(eq. (22a) in ref. 21 ) . 


Composite Methods 

Two composite methods, or methods which rely on one or more second-order 
calculations to achieve their higher order accuracy, were studied. They are 
Richardson extrapolation and the method of deferred corrections applied to a 
second-order method; both are discussed by Keller (ref. 1). The second-order 
method used in the composite method was chosen from among the Keller box scheme 
(KBS) , the conservative Keller box scheme (CKBS) , and the Davis coupled scheme 
(DCS) • The CKBS was chosen because it provides more accurate results than the 
DCS on a coarse grid and is more efficient than the KBS for laminar flows. The 
three second-order methods are compared in detail in the section "Selection of 
Second-Order Method." 

Conservative Keller box scheme with Richardson extrapolation (CKBSRE) .- 
Application of Richardson extrapolation to obtain higher order accuracy using 
the second-order Keller box scheme for the boundary-layer equations has been 
reported by Keller and Cebeci (ref. 22) . Richardson extrapolation is applied 
here to the second-order conservative Keller box scheme (CKBS) , which is the 
second-order method obtained when the higher order terms are deleted from the 
B4S. The extrapolated solution f^ for a grid with I intervals is 
obtained by the following relation: 

fe(l) = f2n-l(2l) + ^[f2n-l (2D - fn(D] (32) 

where solution on a grid with I intervals and fj^(2l) is the 

solution on a grid with 21 intervals. Identical extrapolations are used for 
the dependent variables s^ and Vj^ as well as for the transformed displace- 
ment thickness 6*. Keller (ref. 1) has shown that for second-order methods, 
a single Richardson extrapolation should result in a solution with accuracy 
0 ( 1 - 4 ) . 

Method of deferred corrections (MDC) .- The method of deferred corrections 
is combined with the CKBS by applying equation (11) in the form 

9n - 9n-l “ “(9n + Sn-l > = (33a) 


where 


R 


n 



" V 

9n-l) 


(33b) 


13 



To obtain a solution with the MDC, a converged second-order solution is 
first obtained with = 0. Hie solution thus obtained is used to evaluate 
Rn using the functions for g" given in equations (21) to (23). A second 
converged solution is then obtained for equation (33a) with R^ fixed at the 
value obtained from the first calculation. 

This formulation of the method of deferred corrections requires only one 
deferred correction to make the truncation error formally fourth order which is 
consistent with the B4S. When applying the method of deferred corrections, 
there are basically two approaches in approximating the higher order correc- 
tions; both are discussed by Keller (ref. 1). 

The usual approach is to approximate the derivatives gj^ and gJi-i by 
three-point central differences; thus the values of 9n+l / 9n» 9n-l » 

9n-2 3*^® required to evaluate R^* For uniform grids, the truncation error 

of this approach is formally fourth order for one deferred correction, whereas 
for a variable grid, the formal truncation error is increased to third order. 
This approach also leads to difficulties at the boundaries, where gN+l and 
g _1 are required to compute Rjj and R 2 . Pereyra (ref. 23) avoided this 
difficulty by using one-sided differences at the boundaries for the correction 
terms. Keller discusses a simpler, and probably superior, procedure for treat- 
ing the corrections near the boundaries; Extend the second-order numerical 
solution beyond the boundaries to generate values for 9N+1 ' 9 n+2/ • • • 3*^^^ 

g_ 1 , g_2f • • ♦, as needed for central differences near or at the boundaries. 

An alternative approach for approximating the higher order correction 
terms, and the one adopted here, is to differentiate the governing equations 
to obtain gj^ as a function of the dependent variables f^, s^, and v^. 

While this procedure loses its attractiveness when more than one deferred cor- 
rection is made (i.e., accuracy greater than O(h^) is desired), it appears to 
be more efficient for one correction because no modifications are needed at the 
boundaries and it also has the advantage of remaining formally fourth-order 
accurate for variable grids. 


CALCULATION OF DISPLACEMENT THICKNESS 
The transformed displacement thickness, 

6* = C (1 - f) dn (34) 

'Jq 

is computed for the higher order intrinsic methods and the MDC by using equa- 
tion (11) where g' = 1 - f and thus g" = -s. Equation (11) is integrated 

from the solid boundary, where g-j =0, to the outer boundary, where gjj = 6*. 
The displacement thickness for the CKBSRE is obtained in a similar manner with 
the higher order terms omitted. 
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RESULTS AND DISCUSSION 


The author carefully programmed each method with the same efficiencies to 
obtain a fair comparison with respect to computer run times. The method used 
to solve the 3x3 block- tridiagonal finite-difference equations was coded with 
zero elements eliminated. Several model problems were solved with each method. 
Comparisons between the second-order methods and between the fourth-order meth- 
ods were made for laminar cases with 3 = -0.1, 0, and 1, and comparison between 
the fourth-order methods was made for similar and nonsimilar turbulent cases 
with 3=0 and 0.5. The same initial profiles were used with each method. The 

initial profiles are given by f ^ = 1 - e Sj^ = a(l - fn) r and v^^ = -n^ 

where a = 0.47, 0.47, and 1.23 corresponding to the three values of 3 for 
the laminar cases and a = 2.75 and 3.42 corresponding to the two values of 3 
for the turbulent cases. 

Since a closed-form solution to the governing equations is not available, 
the "exact" solution was taken to be the B4S solution with 640 intervals 
(641 grid points) across the boundary layer, unless otherwise noted. 


Selection of Second-Order Method 

In selecting the second-order method to use with Richardson extrapolation 
and the method of deferred corrections, three second-order methods were studied. 
They are the Keller box scheme (KBS), the conservative Keller box scheme (CKBS), 
and the Davis coupled scheme (DCS). The CKBS is the second-order method 
obtained when the higher order terms are deleted from the B4S. It differs from 
the KBS (ref. 13) mainly in the form of the momentum finite-difference equation 
which is written in conservation form while the KBS version is obtained from 
equation (1 ) which is only partially conservative. The dependent variables for 
both methods are f, s = f*, and v. The DCS (ref. 12) uses a three-point 
variable-grid scheme reported by Blottner (ref. 24) with the linearized finite- 
difference equations for f and v being solved with a modified tridiagonal 
scheme. All three methods solve the momentum and continuity equations coupled 
and are linearized via Newton's method. 

As noted by Keller and Cebeci (ref. 22), the purpose of higher order meth- 
ods is not to acquire more significant digits in the solution but rather to get 
reasonable answers with relatively few grid points. With this purpose in mind, 
the second-order method selected was the one which gave the best results with 
the fewest number of grid points for model problems having favorable as well as 
adverse pressure gradients. 

The second-order methods are compared for laminar-flow problems only. The 
governing equations for the model problems reduce to ordinary differential equa- 
tions and are obtained by deleting the 3 ( )/3^ terms from the governing equa- 
tions and setting A = 1 . The problems were solved with the following numbers 
of intervals across the boundary layer: 10, 16, 20, 32, 40, 64, 80, 128, 160, 

and 320. 
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The grid size for all calculations varied across the boundary layer. The 
equation used to generate the grid is given by 


(1 + 


(35a) 


where 


n - 1 



and C] and a are arbitrary constants chosen here so that for N = 11, 

Am =0.05, and the remaining points are distributed in a reasonable manner. 

For the laminar calculations where Tijg = 24.2538, c-j = -0.4, and a == 8.26. 

The constant cq is obtained by applying equation (35a) at n = N. This type 
of grid stretching was selected because all coarse grid points are also points 
on the exact- solution grid (N * 641); thus, no interpolation between node points 
was required to compute the norms of the f- and s-profiles which were examined 
as a measure of the numerical error present in the finite-difference solutions. 
The f-norm is defined as 


f - f, 


exact 




1 N-1 


N - 2 n=2 


^n “ ^n, exact 


(36) 


with the s-norm being defined in a similar manner. The numerical errors in the 
wall shear and in the transformed boundary-layer displacement thickness 

E-. were also considered in determining the relative accuracy of the methods. 

6 


Figure 1 shows the percentage error in the wall shear versus the number of 
intervals across the boundary layer, I = N - 1, for the CKBS, KBS, and DCS for 
laminar, flat-plate (3 = 0) flow. (The slopes indicated in the figures (2 to 1 
in fig. 1) are the slopes of the lines on the logarithmic grid.) The figure 
shows that as the grid is refined, the CKBS is the most accurate when accuracy 
is measured solely by the percentage error in the wall shear. As the grid is 
refined, the numerical errors for laminar flow have the following form: 


E = CI"^ as I ^ CO 


( 37 ) 


where co = 2 for second-order methods and O) = 4 for fourth-order methods, 

with E being the f-norm, s-norm, %, or E-^ . The constants C were com- 

6 

puted for both 3=1 and 3=0. Typical results show that the ratios of the 
number of intervals to achieve an equivalent accuracy in the wall shear, 
Ikbs/^DCS ^CKBS/^DCS' were 2.3 and 0.5, respectively, for 3=0 and 1.4 

and 0.8, respectively, for 3=1. 
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To compute the relative efficiency of the methods, it was observed that, 
as the grid was refined, the central-processor-unit time in seconds T could 
be approximated as 

T = iCI (38) 

where C is a constant and i is the number of iterations required for con- 
vergence. The efficiency of method A relative to method B is given by 


WTb = (Ca/Cb) (Ia/Ib) = (Ca/Cb) (Ca/Cb)^/^^ (39) 

where Ia/^B ratio of the number of intervals to achieve an equivalent 

accuracy in one of the errors. Because of the short run times, on the Control 
Data CYBER 175 computer, special care was taken in determining the constants C 
with the values used being the average over 6400/1 calculations. Thus, for 
10 intervals, the computer times were averaged over 640 calculations. 

For comparison of the three second-order methods, the convergence criterion 
required that 


|Af| 


max 


< 10-T2 


(40a) 


where 


I I max “ ^n ^n 


(40b) 


Since all the second-order methods converged quadra tically for laminar flows 
and thus required an equal number of iterations for the cases where time effi- 
ciency was computed, it was not felt that this restrictive convergence crite- 
rion penalized any method. Typical results showed that to compute an equivalent 
accuracy in the wall shear, the relative times Tj^bs/'^DCS ^CKBS/^DCS 

were 2.8 and 0.7, respectively, for 3=0. Similar results were obtained for 
3=1. However, the relative efficiency of the methods varied with the pres- 
sure gradient parameter 3 and also with the different numerical errors. 

Note from equation (39) that the relative time efficiency depends both on 
the relative time to solve the difference equations, Cj^bs/^dCS' rela- 

tive truncation error (C^bs/^DCS) ^ number of intervals becomes 
small, figure 1 shows the truncation error to deviate from its asymptotic 
behavior with the DCS showing the greatest deviation. This is borne out in 
table II where the f-profile obtained on a coarse grid (N = 11) is compared 
with the exact solution. 
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TABLE II.- COMPARISON OF SECOND-ORDER f-PROFILES WITH 


EXACT SOLUTION FOR 3=0 
[Laminar flow; N = 11] 




f-profiles calculated with - 


"0 

Exact 

CKBS 

KBS 

DCS 

0 

.49975 X 10-1 

0 

.23468 X lo-l 

0 

.23426 X 10^1 

0 

.24727 X lo-'^ 

0 

.24432 X 10-1 

.14205 X loO 

.66705 X 10“1 

.66585 X 10-1 

.70282 X io-'> 

.69448 X 70-1 

.30761 X iqO 

.14437 X iqO 

.14411 X loO 

.15208 X -[qO 

.15033 X loO 

.60231 X 10° 

.28164 X 10^ 

.28113 X 10° 

.29632 X 10^ 

.29355 X loO 

,11266 X 10 

.51474 X ^Q0 

.51364 X loO 

.53817 X ^oO 

.53930 X IqO 

.20651 X 10 1 

.83286 X 10® 

.83263 X IflO 

.85435 X loO 

.89294 X loO 

,37657 X 10 

.99552 X 10° 

.10140 X 10 

.10156 X 10 

.10820 X 10 

.69004 X 10 

.10000 X 10 

.99643 X loO 

.99600 X loO 

.99698 X 10° 

.12809 X iq2 

.10000 X 10 

.10032 X 10 

.10036 X 10 

.10687 X 10 

.24254 X io2 

.10000 X 10 

.10000 X 10 

.10000 X 10 

.10000 X 10 

Ej, percent . . . 

E-^y percent . . . 

6 


0.179 X IqO 
.167 X 10 

0.537 X 10 
.334 X 10 

0.411 X 10 

.674 X lo2 


For a coarse grid^ the CKBS is slightly more accurate than the KBS. Note the 
very large error (67 percent) in the displacement thickness for the DCS. This 
error appears to result from an approximate 8-percent overshoot in the f-profile 
as compared with 1.5 percent for the CKBS and the KBS. 

The solutions obtained with the three-point DCS on a coarse grid (N = 11 ) 
showed a sensitivity to the pressure gradient parameter B that was not 
observed with either the CKBS or the KBS. For B = 1, the f-profile for all 
three methods showed negligible overshoot with errors in the wall shear and 
displacement thickness that were comparable. However, for B = -0.1, the 
errors in the wall shear and displacement thickness for the DCS were 45 and 
185 percent, respectively, apparently because of a 21 -percent overshoot in the 
f-profile. For this same case, the CKBS and the KBS showed negligible overshoot 
with the error in the wall shear and displacement thickness being 1 percent and 
5 percent, respectively, for the CKBS. 

On the basis of these laminar results on a coarse grid, the CKBS and the 
KBS were judged preferable to the DCS for use with Richardson extrapolation and 
the method of deferred corrections when solutions with only a few grid points 
are desired. Since the CKBS was found to be more efficient than the KBS for 
laminar flows, the CKBS was chosen to use with the higher order composite 
me thods. 


Higher Order Methods Applied to Laminar Flows 

Figure 2 shows the f-norm versus the number of intervals across the bound- 
ary layer for the CKBSRE, MFE, MDC, and RKSl compared with the B4S values for 
laminar, stagnation-point (B = 1 ) flow. 
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Solutions for model problems using the CKBSRE were first obtained on a 
grid with 2l intervals and then on a grid with I intervals* The initial pro- 
files for the I-interval solution were taken from the solution with 2l inter- 
vals* The error norms for the CKBSRE are shown plotted versus the I inter- 
vals* The starting profiles for the MDC calculation where ^ 0 were 

those obtained from the = 0 solution* The exact solution (I = 640) used 
to compute the f-norm for the RKSl method was obtained with that method, since 
it required solution of a slightly different boundary problem (sjg = 0) * 

For laminar flow, figure 2 shows all the higher order calculations with a 
variable grid to be fourth-order accurate as the grid size is refined. The 
f-norm for the CKBS is shown for comparison. As the grid is refined, the error 
norms approach the form given in equation (37) with U) = 4. Numerical values 
for the constants C for this case and the laminar, flat-plate (3=0) case 
were computed* 

To determine the relative time efficiency of the higher order methods to 
achieve equivalent accuracy, the central-processor-unit time in seconds was 
expressed as in equation (38) , The time efficiency was obtained in the same 
manner as for the second-order methods (eq. (39)), 

In a previous paper (ref* 25), the convergence criterion of 
lAfUax^ 10-12 was also applied to the higher order methods when computing 
the relative time efficiency. Further studies have shown that this restrictive 
criterion does in fact penalize the RKSl, MDC, and CKBSRE relative to the B4S. 

To obtain a fair comparison, it is necessary only to require that the con- 
vergence error be less than the truncation error. This criterion was applied 
by noting the minimum number of iterations required for the error to plot as 
shown for I S 128 in figure 2. Results for 6=0 are shown in table III. 


TABLE III.- MINIMUM NUMBER OF ITERATIONS REQUIRED 
FOR CONVERGENCE FOR 6=0 
[Laminar flow; I ^ 128] 


Method 

Number of iterations 

B4S 

5 

RKSl 

4 

MFE 

5 

MDC 

^4 + 1 

CKBSRE 

1 

1^5 + 2 


®The MDC required four iterations to obtain 
the Rjj = 0 solution and one iteration with the 
correction term added. 

^The CKBSRE required five iterations on the 
fine grid and two iterations on the coarse grid. 
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TABLE IV.- RELATIVE EFFICIENCY OF HIGHER ORDER METHODS FOR 6=0 


[Laminar flow] 


Relative times 

For equivalent accuracy 

in - 

Wall shear 

6* 

f-norm 

^MDC/^B4S 

2.3 

0.9 

1.4 

'^CKBSRE/Tb4S 

2.3 

1.1 

2.0 

TmfeAb4S 

2.7 

1.4 

1 .9 

%KS1 Ab4S 

3.6 

1.8 

2.0 


Table IV, which shows the relative time efficiency for this case indicates that 
the B4S is the most efficient higher order method for laminar flows as the grid 
is refined. Similar results were obtained for 6=1. 

As the number of grid points is reduced, the truncation errors deviate 
from their asymptotic behavior. Figure 2 shows the MDC and the RKSl to deviate 
the greatest for the coarse grid, N = 11. Tables V to VII show the f-profiles 
obtained on this coarse grid along with the exact solution for laminar flows 
with favorable as well as adverse pressure gradients (6=1^0, and -0.1). 


TABLE V.- COMPARISON OP HIGHER ORDER f~PROFILES WITH EXACT SOLUTION FOR 6 = 1 

[Laminae flow; N = 11] 


f-profiles calculated with - 



Exact 

B4S 

CKBSRE 

MDC 

RKSl 

MFE 

0 



0 



0 



0 



0 



0 



0 



.49975 

X 

10-^ 

.60350 

X 

10"'' 

.60359 

X 

lO-l 

.60350 

X 

lO'"* 

.60375 

X 

10-1 

.60268 

X 

10“1 

.60376 

X 

10-1 

.14205 

X 

IQO 

.16503 

X 

10° 

.16506 

X 

loO 

.16503 

X 

10^ 

,16510 

X 

10° 

.16480 

X 

10° 

.16510 

X 

10^ 

.30761 

X 

lOO 

.33236 

X 

IQO 

.33242 

X 

loO 

.33236 

X 

loO 

.33251 

X 

10^ 

.33186 

X 

10^ 

.33254 

X 

10^ 

.60231 

X 

lOO 

.56784 

X 

IqO 

.56795 

X 

lOO 

.56787 

X 

10^ 

,56813 

X 

10° 

.56683 

X 

loO 

.56815 

X 

IqO 

.11266 

X 

10 

.82381 

X 

10® 

.82393 

X 

iqO 

.82399 

X 

10^ 

,82400 

X 


.82158 

X 

10° 

.82432 

X 

iqO 

.20651 


10 

.97721 

X 

10® 

.97683 

X 

IQO 

.97725 

X 

iqO 

,97588 

X 

loO 

.96919 

X 

loO 

.97409 

X 

10° 

.37657 

X 

10 

.99989 

X 

10° 

.99904 

X 

IqO 

.99948 

X 

10° 

.10029 

X 

10 

.99194 

X 

10° 

.10016 

X 

10 

.69004 

X 

10 

.10000 

X 

10 

.99972 

X 

loO 

.10003 

X 

10 

,10018 

X 

10 

.99934 

X 

IqO 

.99965 

X 

10^ 

.1 2809 

X 

102 

.10000 

X 

10 

.99984 

X 

IQO 

.99985 

X 

loO 

,98816 

X 

loO 

.92008 

X 

10° 

.10002 

X 

10 

.24254 

X 

102 

.10000 

X 

10 

.10000 

X 

10 

,10000 

X 

10 

.10000 

X 

10 1 

.10000 

X 

10 

.10000 

X 

10 

Ex, percent 

. . . 




0.153 X 

1C 

rT 

0.612 ^ 

10 

1-4 

0.405 X 

1C 

rl 1 

0,132 X 

lOO 

0.443 X 

1C 

1-1 

E- . , percent - . . 
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.177 X 

10° 

.104 X 

10° 

1 

.334 X 

1Q2 

.129 X 

lo2 

.166 X 

loO 
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TABLE VI. - COMPARISON OF HIGHER ORDER f-PROFILES WITH EXACT SOLUTION FOR 3=0 

[Laminar flow; N = 111 


n 

0 

.49975 X 10-1 

Exact 

0 

.23468 X 10-1 

B4S 

0 

.23480 X 10-1 

f-profiles cal< 
CKBSRE 

0 

.23456 X 10-1 

:ulated with - 
MbC 

0 

.19602 X 10-1 

RKSl 

0 

-.62917 X 10-3 

MFE 

0 

.23478 X 10-1 

.14205 X ^oO 

.66705 X 10-1 

.66738 X 10-1 

.66671 X 10-1 

.55715 X 10-1 

-.17884 X 10-2 

.66724 X 10-1 

.30761 X iqO 

.14437 X loO 

.14444 X loO 

.14430 X IqO 

.12060 X loO 

-.38729 X 10-2 

.14445 X 10^ 

.60231 X iqO 

.28164 X loO 

.28178 X IqO 

.28149 X loO 

.23541 X loO 

-.75839 X 10-2 

.28166 X loO 

.11266 X 10 

.51474 X 10° 

.5T498 X loO 

o 

o 

X 

in 

.43191 X loO 

-.14194 X 10-1 

.51517 X loO 

.20651 X 10 

.83286 X loO 

.83293 X loO 

.83278 X loO 

.71134 X 10° 

-.26121 X 10-1 

.83429 X ^oO 

.37657 X 10 

.99552 X loO 

.99616 X loO 

.99536 X loO 

.87000 X loO 

-.48672 X 10-1 

.99275 X loO 

.69004 X 10 

.10000 X 10 

.99947 X loO 

.10012 X 10 

.91739 X loO 

-.11413 X loO 

.10003 X 10 

.12809 X 102 

.10000 X 10 

.99983 X loO 

.99893 X loO 

.83792 X loO 

-.34076 X loO 

.99871 X loO 

.24254 X lo2 

.10000 X 10 

.10000 X 10 

.10000 X 10 

.10000 X 10 



.10000 X 10 

.10000 X 10 

Ej, percent . - . 


0.500 X 10-1 

0.496 X 10-1 

0.165 X io2 


0.483 X 10-1 

percent . . . 


.345 X loO 

.714 X loO 

___ . 

.231 X lo3 

.. 


o 

o 

X 

CM 

00 

CM 


TABLE VII.- COMPARISON OF HIGHER ORDER f-PROFILES WITH EXACT SOLUTION FOR 3 = -0.1 

[Laminar flow; N = ll] 











f-profiles 

calculated with - 










Exact 


B4S 


CKBSRE 

MDC 

RKSl 

MFE 

0 



0 



0 



0 



0 



0 



0 



.49975 

X 

10-1 

.16080 

X 

10-1 

.16092 

X 

10-1 

.16083 

X 

10-1 

.14291 

X 

10-1 

-.13682 

X 

10“2 

.16002 

X 

10-1 

.14205 

X 

loO 

.46360 

X 

10-1 

.46394 

X 

10-1 

.46368 

X 

10-1 

.41276 

X 

10-1 

-.32351 

X 

10-2 

.46195 

X 

10-1 

.30761 

X 

loO 

.10289 

X 

loO 

.10297 

X 

loO 

.10291 

X 

10° 

.91894 

X 

10-1 

-.44593 

X 

10-2 

.10228 

X 

lOO 

.60231 

X 

loO 

.20967 

X 

loO 

.20981 

X 

10° 

.20970 

X 

10° 

.18827 

X 

lOO 

.14425 

X 

10-3 

.20937 

X 

10° 

.11266 

X 

10 

.41286 

X 

loO 

.41312 

X 

lOO 

.41284 

X 

loO 

.37455 

X 

loO 

.29807 

X 

10-1 

.40926 

X 

lOO 

.20651 

X 

10 

.75159 

X 

loO 

.75173 

X 

10° 

.75152 

X 

o 

o 

.69731 

X 

loO 

.15007 

X 

10° 

.75645 

X 

loO 

.37657 

X 

10 

.98964 

X 

loO 

.99005 

X 

loO 

.99217 

X 

lOO 

.94324 

X 

io?J 

.50546 

X 

10° 

.98016 

X 

loO 

.69004 

X 

10 

.10000 

X 

10 

.99900 

X 

lOO 

.99980 

X 

loO 

.97997 

X 

loO 

.79147 

X 

loO 

.10000 

X 

10 

.12809 

X 

102 

.10000 

X 

10 

.99971 

X 

o 

o 

.99964 

X 

lOO 

.96007 

X 

loO 

.82942 

X 

loO 

.99682 

X 

lOO 

.24254 

X 

102 

.10000 

X 

10 

.10000 

X 

10 

.10000 

X 

10 

.10000 

X 

10 

.10000 

X 

10 

.10000 

X 

10 

Et# percent 

. . . 




0.743 X 

10-1 

0.172 X 

10 

-1 

0.112 X 

102 

0.109 X 

lo3 

0.563 X 

10° 

E-. , percent . . . 
6 




.384 X 

O 1 

o 

.864 X 

10 

-1 

.538 X 

102 

.322 X 

103 

.129 X 

loO 
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when both displacanent thickness and wall shear are considered, table V 
(B = 1) shows the accuracy of the CKBSRE, MFE, and B4S to be about the same, 
while the MDC and RKSl both show significant errors in the displacement thick- 
ness. Note that the MDC and the RKSl show large oscillations in the f-profile. 
Table VI (B = 0) shows similar results for the MDC, while the RKSl method gives 
a nonphysical solution. Results for an adverse pressure gradient (B = -0.1) are 
shown in table VII where behavior similar to the other cases (B = 1 and B = 0) 
is observed. 


^^plication of Higher Order Methods to Turbulent Flow 

Similar calculations .- Figures 3 and 4 show the percentage error in the 
wall shear versus the number of intervals across the boundary layer for calcu- 
lations of a model turbulent problem. The governing equations for these cases 
are obtained by deleting the 3 ( )/3^ terms from the governing equations. The 
Reynolds number R^ = 1.88 x 10®, ^ = 1, and Ug = 1 , with the outer boundary 
located at Hu = 60. The constants in the grid-stretch equation (35) are 
c-| =0.05 and a = -109 and were selected so that for N = 11, Arj-j =0.05 
and the other points were stretched to hu =60 in a reasonable manner. The 

scale length 7 = 24.5 for B = 0 and 7 = 15.8 for B =0.5. The con- 
vergence criterion for these calculations is given by equation (40a). To 
achieve this convergence, the MDC required multiple iterations with the cor- 
rection term added. 

For B = 0, figure 3 shows the B4S, MFE, and the CKBSRE to be significantly 
more accurate for the same grid than the MDC and the RKSl , both of which are 
less accurate than the second-order CKBS for I < 128. The converged solutions 
shown for the RKSl were obtained with underrelaxation (relaxation factor 
r = 0.5). An exact solution (i.e., using 641 grid points) with RKSl could not 
be obtained, and thus, the numerical errors were computed using the exact solu- 
tion obtained with the B4S. 

For B = 0.5, figure 4 shows the B4S, MFE, and CKBSRE again to be more 
accurate than the MDC. Converged solutions using RKSl could not be obtained 
for the grid sizes shown even with underrelaxation factors r as small as 0.5. 
However, an exact solution with RKSl was obtained (r = 0.7). The exact wall 
shear values computed with B4S and RKSl agreed to six decimal places. 

Because of the roughness of the error curves shown in figures 3 and 4, it 
was not feasible to make time efficiency estimates similar to those for laminar 
flows. The reason for this error scatter in the turbulent cases is not known. 
However, similar behavior is evident in a plot given by Blottner (ref. 24) . 

The derivative f" was observed to be smooth for this case and all other cal- 
culations with the B4S. 

Nonsimilar calculations .- To determine the relative efficiency of higher 
order methods for nonsimilar turbulent calculations, the B4S, CKBSRE, MFE, and 
CKBS were applied to compute the incompressible turbulent boundary layer on a 
flat plate. The C-grid is given by ^ = 0, 0.00086, 0.0043, 0.0086, 0.043, 

0.17, 0.26, 0.34, 0.51, 0.68, 0.86, 1.03, V.28, 1.50, 1.71, and 1.88. The 
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reference Reynolds number 1 x 10^ giving a local Reynolds number 

Rjj = 1.88 X 10^ at C = 1.88. The ri-grid was given by equation (35) with 
tIn = 24.2538. The constants c^i and a in the grid-stretch equation were 
-0.4 and 8.26, respectively. 

Second-order solutions to this problem have been reported by Keller and 
Cebeci (ref. 22), who also presented higher order results using Richardson 
extrapolation (An 0) , and Blottner (ref. 24). Both of these calculations 
were made using a two-layer eddy viscosity formulation and a grid-point distri- 
bution such that Arij^/Arin-i - Constant. The grid-point distribution used here 
closely resembles that used in those previous calculations, that is, Af]] = 0.05 
for N = 11, while a single-layer eddy viscosity model is used for the present 
solutions. 

Figure 5 shows the percentage error in the wall shear versus the number of 
intervals across the boundary layer at C = 1.88. At 5 = 0.00086, the Crank- 
Nicolson scheme was used to approximate the 3 ( )/3C terms, while at all other 

stations the three-point backward difference scheme was used. The value of l 
was 15.8. Transition was at K - 0.00086. 

The exact solution was obtained with 640 intervals across the boundary 
layer with the C-grid given previously. Richardson extrapolation for A^ 0 
was not made and, as such, the term "exact" means only as the Ari grid is 
refined. The convergence criterion used to obtain the results shown in fig- 
ure 5 is given by equation (40a) . The only exceptions are for the MPE and B4S 
where fully converged solutions at the downstream stations could not be reached 
for I ^ 20 for the MPE and 1*10 for B4S. 

For laminar calculations, the error norms are smooth and vary according to 
equation (37) for o) = 4. Thus, comparisons for relative time efficiency of 
the higher order methods can be computed using equation (39) and the results, 
shown in table IV, apply over a wide range of intervals (I < 128). However, 
for turbulent calculations the error norms are not well behaved and, thus the 
relative efficiency cannot be computed over a range of intervals as for the 
laminar calculations. 

For nonsimilar turbulent calculations, the efficiencies of the CKBSRE, 

MPE, and CKBS relative to the B4S were computed for a specific accuracy at a 
specific point, that is, for a 0.1 -percent numerical error in the wall shear at 
5 = 1.88. The convergence criterion for these comparisons required that the 

absolute change in 6* between consecutive iterations be less than 1 percent. 
Shown in table VIII are the results where, to obtain this numerical accuracy, 
1=12 for both the B4S and CKBSRE and I = 20 and 70 for the MPE and CKBS, 
respectively. Both the MPE and B4S converged without trouble to this criterion, 
requiring three iterations per streamwise station. For this example, the MPE 
is a factor of 2 less efficient than the B4S. Note that the second-order CKBS 
is less efficient than the B4S by a factor of 4.8. 
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TABLE VIII.- RELATIVE EFFICIENCY OF HIGHER ORDER METHODS FOR 
NONSIMILAR TURBULENT CALCULATIONS 


Method 

Number of intervals 

Relative efficiency^ 

B4S 

12 

1 .0 

MFE 

20 

2.0 

CKBSRE 

12 and 24 

2.5 

CKBS 

70 

4.8 


^Relative time efficiency with respect to B4S to achieve 
0.1 -percent numerical error in wall shear at ^ = 1.88, 

Ryr =1.88 X 10^. 


CONCLUSIONS 

A fourth-order box method for calculating numerical solutions to parabolic, 
partial-differential equations in two variables or ordinary differential equa- 
tions has been presented. The method, which is the natural extension of the 
second-order box scheme to fourth order, has been demonstrated with application 
to the incompressible, laminar and turbulent, boundary-layer equations. This 
method has been compared with other two-point and three-point higher order 
methods on variable grids for laminar and turbulent, incompressible flows. On 
the basis of these numerical results the following conclusions were reached: 

(1) For laminar calculations, the present fourth-order box scheme (B4S) , 
the conservative Keller box scheme with Richardson extrapolation (CKBSRE) , the 
method of deferred corrections (MDC) , the modified finite-element method (MFE) , 
and a three-point spline method are fourth-order accurate with variable grids. 

(2) For laminar calculations on coarse grids, the B4S, MFE, and CKBSRE 
give greater accuracy than the MDC and the three-point spline method. 

(3) For laminar calculations, the B4S is the most efficient higher order 
me thod. 

(4) For turbulent calculations, the B4S, MFE, and CKBSRE are significantly 
more accurate than the MDC or the three-point nonconservative spline method. 

(5) For a nonsimilar turbulent calculation, the B4S showed a factor of 4.8 
reduction in computer time to achieve the same accuracy as the second-order 
conservative Keller box scheme, a factor of 2 reduction compared with the MFE, 
and a factor of 2.5 reduction compared with the CKBSRE. 
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(6) For practical two-dimensional, incompressible, turbulent boundary- 
layer calculations (i.e., numerical errors of approximately 0.1 percent), the 
B4S is the most efficient higher order method. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
September 14, 1978 



APPENDIX A 


COEFFICIENTS FOR B4S 

In this appendix, it is shewn how to obtain the elements a^j and b^j 
of the matrices A^ and and the components p^f <3nr ^n the vec- 
tor which appear in equations (24) . Equations (24) are the matrix equa- 

tions which result when equations (15) to (23) are substituted into the finite- 
difference expression (eq. (11)) and these finite-difference equations are 
linearized . 

The 9( )/3? terms in equations (18) and (20) to (23) are written in 
finite-difference form as 


3 ( ) 


m- ( 1 -6 ) , n 


= Cl ( 


)m,n ” °2( )m-l,n °3 ( )m-2,n 


(Al) 


where m,n are the grid indices in the ^ ,r\ coordinates. The 3( )/3C terms 
are written so that they may be approximated either by the Crank-Nicolson scheme 
(6 = 1/2) where 


1 


Cl = C2 = 


S3 = 0 


^^m-1 

or by a three-point backward-difference scheme (6=1) where 

3 


(A2) 


Cl = 


(1 + Y)^^m-2 


(A3) 


2 


(A4) 


C3 


2y - 1 

(1 + Y)^^m-2 


with 


A?m-1 

ACm-2 


^^m-1 “ ^m “ ^m-1 


(A5) 


(A6) 


(A7) 
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The similarity form of the difference equations is obtained by setting 
6=1 with Cl = C 2 = S 3 = 0. 


Momentum Finite-Difference Equation 
Equation (15) is written as 


9n ~ 0 (®'m,n®m,n 


'^mjn^mfn)^ + (1 - ® ) (^-m-l ,n®m-l ,n ” ^m -1 ,n^m-l ,n) 


(A 8 ) 


where i denotes the present iteration. Nonlinear products such as 
(Vjn^n^m,n)^ are linearized by Newton's method which can be shown to result in 


' -in ,1 




The linearized form of equation (A 8 ) becomes 


(A9) 


9n = 0 


( i /n j. / « \i-l 

\2*m,n ” V®m»n ~ ^^m,n " ‘J^,n ~ %,n'm,n “ ^^m^n^mjn ■*" ('^m,nfm,nl 


+ (1 - 6 ) (S,jj_i ^nSm -1 ,n “ ^m -1 ,n^m-l ,n) 


The values of ^ and 3 are defined as 


(AlO) 




3 = 63n, + (1 - 6)3n,-i 


(All) 

(A12) 


Equation (18) for g/, becomes 


9n “ °1 ^^m,n) ~ ^ 2 ^m-l ,n "*■ ®3^m-2,n^ + (1 + 


+ (1 - 0 ) fm-i ,n 


- 3 


(A13) 


Now define 


5i = 2Cci + (1 + 3)6 


(A14) 
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«2 = 

-2^2 + (T + 3) (T - 6) 

(ATS) 

Sb = 

25c3 

(AT 6) 


The nonlinear expressions for and become 

9n ~ ,n ®3%-2,n “ ® (A17) 

and 


Jn = 2[ 


9n ~ (fm,nSm,n) ^ "*■ “ 2 ^ 111-1 ,n®m-l ,n "*■ ®3fm-2,n®i 


m-2,n_ 


(A18) 


After linearizing equations (AT 7) and (AT 8) via Newton’s method, the lin- 
earized expressions for g„, g^, and g^ are substituted into equation (TT), 
and the resulting equation is rewritten in the following form: 


3TTS^,n-T + bns^^n + aT2fm,n-T + bi2fm,n + aT3Vm,n-T + bi3Vm,n = Pn 
where 

3tt = “ ~ ®T^m,n-T 

bTT = ©(sS-mTi - l) 

o i-1 up: ^i-1 i-1 

3T2 = ^Vjn^n-T " haifj„^n-T “ — <^TSin,n-T 

D 

5T2 = “ ^T^m,n ■*■ ~ ^T%,n 

D 

3T3 = ®fin,n-T 
*^T3 ~ “®^m,n 


(AT 9) 

(A20) 

(A2T) 

(A22) 

(A23) 

(A24) 

(A25) 
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Pn "" 


Pn = 


and 

loop 


n " ” ^nirn^rn “ ^nif n-1 J 

-Si + (fm,n-l)^’^] + «! ^[ (fm,nSm,n) 

- ^^m,n-l ®m, n-1 ) ^ ^ J ■*■ Pn 


(A26) 


“(1 - 6 ) ,nSm-l ,n ” ^n^m-1 ,n) “ (^m-1 ,n-l Sm-1 ,n-l 

- vm-i ,n-lfm-l ,n-l)] +~[s2(fm-l,n + fm-l,n-l) + S3^f^_2^n + fm-2,n-l) “ 2 bJ 

h^f- 

^ 1^2 (fm-1 ,n®m-l ,n ^m-1 ,n-l %-l ,n-l ) 


^3 (^m-2,n®m~2 ,n " ^m-2,n-l %-2,n-l ) J 


(A27) 


(~) denotes quantities that are evaluated outside of the main iteration 

Shear Finite-Difference Equation 

For the shear equation, gj^ (eq. (16)) and (eq. (19)) are written as 


^n “ ®^m,n ®)^m-l,n 


(A28) 


^n “ ®%,n + " ®)%-l,n 


(A29) 
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To obtain a linearized expression for equation (22) is rewritten in the 

following form: 


^®(^m,n + (1 - ” ®^^m,n®mrn)^ + (1 - ,n^-^ ,n^ 


+ 3 


+ (1 - " 1 


and 


+ K 


Cl 


f 2 ~ 2 ~ 2 

\^m,n/ ~ C2fm-l,n C3fm-2,n 


” ®*in,n ®m,n 


i-1 
®m, n 


Now define 


~ (1 “ ® ) *in-l ,nl ,nl ®m-l ,n (A30) 


«4 = ?c-j + 03 


(A31) 


«5 = —^02 + (1 “ 0)3 


(A32) 


“6 = ^C3 


(A33) 


~ ® ) ^^m-1 ,n®m-l ,n “ ,n| ,n| ,n^ "*■ ®5fm-l,n "*■ ®6^in-2,n “ ^ 


(A34) 


30 


APPENDIX A 


Equation (A30) can be rewritten as 



- 1 ) 





(A35) 


After linearizing equation (A35) via Newton's method, the following linearized 
expression is obtained for gj^t 


9n = 



^m,n®m,n ~ (''m ,n%,n) 


+ a> 


2f f ^ 




where 


Hn = 


e^2S- 


i-l 

m,n 



- 9)(2Jln,.i,n 



(A37) 


Upon substitution of equations (A28) , (A29) , and (A36) into equation (11), the 
shear equation can be rewritten as 


32lSm,n-l + b2is^,n + a22fm,n-l + ^22fm,n + a23Vm,n-l + t.23Vm,n = ^n 


(A38) 


where 



2^m,n-l 



(1-6) ,n-l 



(A39) 
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b2l 


h h2 

- - 6 + e — 

2 12 



(i_0 ) 


i-1 



h2 


.i-1 


322 - “6 +614 — Hn_ifjj,^n-1 
6 


b22 = 6 


~ b2 j_i 

+ ~ Hnfm,n 


323 


h2 

“ ~® TT ®n-l®m,n-l * 


12 


t>23 


h2 

= e — H: 
12 


i-1 

n®m,n 



+ 


26 (j. 


i-1 

m,n-l 



• f‘ 

5m-(l-0),n-lJ 



+ % 


q„ = -(1 - 6) 


^m-1 ,n ~ fm-l,n-l “ „^®m-l,n 


Sm-i ,n-l ) 

— I 


with Kji = ±1, having the same sign as 


i-1 

®m,n* 


(A40) 


(A41) 


(A42) 


(A43) 


(A44) 


1 


(A45) 


(A46) 
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Continuity Finite-Difference Equation 
Equations (17), (20), and (23) can be written as 

9n = + n - 6)Vm-i,n 

®^in,n “ “ ®)fm-l,n 

Define 


9n “ **■ ^3fm--2,n 

" ^2Sm-l ,n ■** ^3%-2^n 


5 7 = “2^c*| - 0 


oJq = 2Cc2 -0-0) 


59 = -2^C3 


Then, substitution of equations (A47) , (A48) , and (A49) into equation (1 
yields the following equation: 

+ “8(fm-l,n + fm-l,n-l) + “9(fm-2,n + fm-2,n-l)] + ~[“?(s^,n ~ %,n-l) 

■'■ ®8(Sm-l,n “ %-l ,n-l ) ®9(Sni-2,n “ Sm-2,n-l)J ® 0 


which can be written as 


a31%,n-l + ^31%,n + a32fm,n-l + t)32fm,n + a33Vm,n-l + b33Vm,n = *^n 


(A47) 

(A48) 

(A49) 

(A50) 

(A51) 

(A52) 


(A53) 

(A54) 
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where 


. h2 
asi = -a, 


(ASS) 


h2 

^31 “ “7 — 

12 


h 

332 = -S? j 


h 



333 = -e 

t>33 = 6 

= ?n 


(AS6) 

(AS7) 

(A58) 

(AS9) 

(A60) 

(A61) 


- e)(vni_T^n “ Vm-l,n-l) 


+ “ “8(fm-l ,n + 


fm-1 ,n-l) 


+ Sg (fni-2,n 


1 h2r 


%- 1 , n- 1 ) 


+ ®9(%-2,n ■ 


%-2,n-l) 


] 


(A62) 
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METHOD OF SOLUTION (B4S) OF FINITE-DIFFERENCE EQUATIONS 

The linear system of finite-difference equations (24) can be solved in the 
relatively simple manner which is presented in this appendix. If equations (24) 
are applied at n = 2 (n = 1 being the surface and n = N locating the outer 
boundary) and the boundary conditions f*| = V] = 0 are applied, three equations 
with four unknowns are obtained. Thus, three of the unknowns may be expressed 
in terms of the fourth unknown. 

Since fjj is known at the outer boundary, it is convenient to solve for 
si , S 2 f and V 2 in terms of f 2 , in anticipation of the general form of the 
recursion relations to be derived subsequently. Thus, at n = 2, equations (24) 
can be rewritten as 


(s2fSi ,V2)’’ = 



(bi2»b22 



where 


bll 

an 

bl3 

b 21 

321 

b23 

b31 

331 

b33 


(Bl) 


(B2) 


Next by applying equations (24) at n = 3 and using equation (Bl ) to 
eliminate S 2 and V 2 as unknowns, the result is again three equations and 
four unknowns, namely, f2# ^3r ^3, and 33. Thus, three of the unknowns may 
again be expressed in terms of the fourth unknown; that is, solve for f 2 ^ V 3 , 

and S3, in terms of f 3 (this order is dictated by the choice at n = 2 which 
was made with the consideration of the outer boundary condition f = 1 ) . If 
this procedure is repeated for n = 4, . . ., N, the following general form is 
obtained for the three unknowns written in terms of the fourth: 


Sn = <3n 


0 ) . _0) 




+ 6r 


(B3) 


_ ( 2 ) ( 2 ) 
^n-1 “ ^n ®n ^n 




(n = 3 , 4 , . , N) 


(B4) 


„ _ fl(3) 4. «(3)^= 


(B5) 


vAiere the superscripts 1,2, and 3 in parentheses identify the coefficients 
associated with the unknowns s, f, and v, respectively. 
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The coefficients ^ f r given 

by 

^^n ^ ^^n ^ ^^n ” ^n'^^^Pn'^Jn^'^n)^ ” t®ll /®21 /®31 “ ^®13/®23»®33)^‘^n-i J 

(B6) 

(B7) 

with 


j = -Dr"^ (bl2rt>22/t>32)^ 


bn ai2 + an ©nil + aisenii bi3 


Dn = 


b2i 322 + a2ienli + 32360-1 b23 


(B8) 


J^31 ®32 + ®31®nM + ®33®n-i ^>33 


Examination of equations (B3) to (B8) reveals that although the solution for the 
dependent variables s^f for and Vj, cannot be computed until the outer bound- 
ary is reached where f jj = 1 , the coefficients <3n^ ^ r <3n^^r *3n^^r ®n^ ^ » 

e^^^, and can be computed if the initial values <32^^/ ®2^ ^ » and 

(3) 

62 are known. These initial values are found by comparing equations (B3) 
and (B5) evaluated at n = 2 with equations (Bl). The initial values are 


f (1) (3)1^ 

1^62 »‘^/®2 J 


C~^ (P2f<32*r2)'*^ 


-C~^ (bi2»b22»b32)^ 


(B9) 


(BIO) 
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vrfiere ^ and w are coefficients to compute s*| from 


s-| = + 0 )f 2 


(BID 


When the outer boundary is reached, equations (B3) to (B5) can be solved 
for Sjj# r and vjj since the outer boundary condition f = 1 can be 

applied. Next sjg-i , fN-2r and vjj-i can be computed since fN-1 is known; 
and likewise all the values of Sj^, and Vj^ proceeding from the outer 

boundary toward the solid boundary can be computed. 
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Intervals 


Figure 1.- Accuracy of wall shear with second-order methods 
Laminar flow; B = 0; s-) = 0.469599988 (exact). 






Intervals 


(a) CKBSRE with B4S. 


Figure 2.- Comparison of higher order methods for laminar flow. 
B = 1; SI = 1.232587657 (exact). 
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(d) MFE with B4S. 
Figure 2.- Concluded. 
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(c) CKBSRE with BHS. 


Figure 3-- Continued. 
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